Analysis of a dissipative model of self-organized criticality with random neighbors 
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We analyze a random neighbor version of the OFC stick-slip model. We find that the mean 
avalanche size is finite as soon as dissipation exists in the bulk but that this size grows exponentially 
fast when dissipation tends to zero. 
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It is an appealing idea that many power laws observed 
in nature arise from an intrinsic trend of a large class 
of extended non-equilibrium systems to evolve toward 
critical points This concept of self-organized criti- 
cality (SOC) has therefore attracted much interest and 
implicit assumptions of the original model have been sub- 
jected to intense scrutiny. Several early SOC explana- 
tions assigned a crucial role to strict bulk conservation 
1^,1). Non-conserving models which show SOC behav- 
ior have since been found j|,[3) but in several cases the 
effect of a small dissipation remains unclear. An interest- 
ing model where dissipation is controlled by a parameter 
a, referred hereafter as the OFC model, has been intro- 
duced in H as a simplified version of previous modelling 
of fault dynamics [Q . Numerical results and supporting 
arguments appear to indicate that the OFC model ex- 
hibits power-law distributed avalanches in the dissipative 
range of a values below the conserving ao (^J^]. A ran- 
dom neighbor version of the OFC model has been studied 
in H and numerical evidence of SOC behavior has sim- 
ilarly been found for u c < a < at®. Our aim here is to 
analyze this simpler version of the OFC model. In con- 
trast to |9) , we find that avalanches are of finite size up 
to the conserving limit a — a® but that their mean size 
grows exponentially fast as a — ► ao. This may explain 
our disagreement with J^] and can perhaps also serve as 
a cautionary note about similar numerical evidence ob- 
tained for local lattice models. 

The model || consists of a set of N sites, to each of 
which is associated an "energy" a;,-. The dynamics alter- 
nates between two phases: 

- the loading phase is supposed to take place on a long 
time scale in the fault dynamics context. In this phase, 
all the below a threshold and increase continu- 
ously and simultaneously with time. This regime lasts 
until one energy reaches the threshold energy which we 
choose equal to one. At this point, an avalanche starts 
and the dynamics enters the avalanche phase. 

- the avalanche phase is thought to be instantaneous on 
the time scale of the loading which can thus be neglected. 
The dynamics is entirely governed by energy transfers be- 
tween different sites. For each xt > 1, K different sites 
j'W are randomly chosen. On each one, the energy is in- 
creased from Xj(i) to Xj(i) + axi and then Xi is set to 0. 



The process is repeated if some of the new energies are 
above one. When all the site energies are smaller than 
one, the avalanche ends and the system returns to the 
loading phase. 

The parameter < a < 1/K controls the dissipation 
during an avalanche. If a = 1/K, energy is conserved 
and the total energy of the system is constant during an 
avalanche. On the contrary, it decreases for a < 1/K. 




FIG. 1. Stationary distribution obtained by simulation 
(circles) and by solving the equations (JsJ) and ([?]) (continu- 
ous line) for K — 4 and a — 0.2. Simulations results are the 
average of 2 10 4 avalanches on a lattice of N = 5 10 3 sites 
and the bin size is 0.01. 



Results of simulations ||] as those reported in Fig.l, 
have shown that the probability distribution Pt(x) of 
site energies tends at large times toward a non-trivial 
stationary distribution P(x). It is, in fact, possible to 
obtain the exact evolution equation obeyed by Pt(x) in 
the limit N — > 00 as we now show. We define the size 
of an avalanche as the number of topplings during its 
evolution (the number of sites where the energy becomes 
greater or equal to the threshold). We suppose that the 
parameter a is small enough so that the mean avalanche 
size has a finite value when the system size tends to in- 
finity. The two regimes of the dynamics contribute to the 
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evolution of Pt and are considered in turn. For definite- 
ness, we fix to unity the growth rate due to the constant 
loading. So, in a time interval At between t and t + At, 
the site energies increase by Ax = At. This gives rise 
to M avalanches with M = NP t (l)Ax (NP t (l) is the 
density of sites which have their energy just below 1). 
The sites which are updated during the course of these 
avalanches belong to three distinct classes: 
-i) the M starting sites of the avalanches the energy of 
which is set to zero, 

-ii) the sites which through energy redistribution have 
their energies first increased above one and then set to 
zero, the total number of which we define to be MA t , 
-iii) those which have their energies increased below one, 
the final energies of which are distributed according to a 
density MB t (x). 

Note that we assume that the system is large enough so 
that the probability that a given site has been updated 
more than once is negligible. At t + At, the probability 
distribution of site energies has become Pt+At(%) with 

NP t+At (x) = NP t {x - Ax) - MS(x - 1) + MB t {x) 

+ M(A t + l)S(x) - KM (A t + l)P t (x). (1) 

The last term on the r.h.s accounts for the fact that the 
sites of classes ii) and iii) are picked at random among 
the iV sites and that their total number is KM (At + 1) 
(since the energy of each site above threshold is redis- 
tributed to K other sites). Taking the limit At — > 0, we 
obtain the evolution equation for Pt(x), < x < 1 : 

d t Pt(x) + d x P t (x) = P t (l) [B t (x) - K(A t + l)P t (x)] (2) 



together with the injection condition at x = 0, 
P t (0)=Pt(l)(It + l) 



(3) 



In order to obtain a closed equation for Pt(x), it re- 
mains to determine the characteristics of the avalanches, 
A t and B t (x), in terms of Pt(x). To this end, we analyze 
the course of an avalanche in a more detailed way. At the 
n-th step of an avalanche, for each site in class ii) which is 
set to 0, K sites are randomly chosen. Those which have 
their energies temporarily increased above one belong to 
class ii) and we denote by a n (x), x > 1, the distribution 
of their energies above threshold. The avalanche ends 
when a n (x) — 0. Similarly, we call b n (x), a < x < 1 
the distribution of sites of class iii) produced at the n-th 
step. This gives therefore for x > 1, 



a n \x 



K 



a n -i(y)P(x - ay) dy 



(4) 



For x < 1, one obtains an equation with the same r.h.s. 
but with b n (x) instead on the l.h.s.. In Eq. (0), the 
integral upper bound has been taken to be 1/(1 — a) 
since a brief study of the series defined by uq — 1 and 
Un+i — l + au n shows that a n (x) is zero if x > 1/(1 — a). 



To compute the evolution of Pt(x), we can restrict our- 
selves to consider B t (x) which is the total density J2 n 
of sites of class iii) averaged over the M avalanches oc- 
curring between t and t + At. We similarly define A t (x) 
as the average over these avalanches of ^2 n a n . From (0), 
A t (x) is determined from Pt(x) as the solution of the 
linear equation for x > 1 : 



A t (x) = K 



A t (y)P t (x - ay) dy + P t (x - a) 



(5) 



For x < 1, the l.h.s. is replaced by B t (x), 



B t (x)=K 



A t (y)P t (x - ay) dy + P t (x - a) 



(6) 



This gives B t in terms of A t (x) and Pt(x). Since A t = 
/ dx A t (x), the evolution of Pt(x) is determined by solv- 
ing Eq. (g) together with (||) using the expression @ for 
B t (x) . Specializing to the steady state time independent 
functions, we finally obtain for < x < 1, 

P'(x) | _ J 1 T ^A(y)P(x-ay)dy + P(x-a) (?) 



KP(0) 



A + l 



where the steady state distribution A(x) is determined 
from P(x) by Eq. (|^) (with the time indices dropped). 
At this stage, several simple remarks can be made. The 
size of an avalanche is the total number of sites in class ii) 
plus the starting site so the mean avalanche size s = 1+A. 
It is useful to note that the injection condition (^) is di- 
rectly obtained by integrating Eq. (|?]) from x = to x = 1 

and by using J \P = 1 and J B = K(A+ 1) - A (the last 
equality follows from the avalanche rule as noted above 
but it can also be derived by adding Eq. (||) and (||) and 
integrating over x). Eq. (||) gives the alternative expres- 
sion of s as s = P(0)/ P(l). There are infinite avalanches 
with non-zero probability only if P(0) is infinite or P(l) 
is zero. Actually, we find below that both are true. Large 
avalanches lead to large P(0) but also deplete the distri- 
bution of sites energies away from a small number of given 
energies with a depletion length proportional to 1/P(0) . 
This leads P(l) to decrease exponentially fast when P(0) 



increases. 



We now turn to the solution of Eq. (H) and (Q). First, 
one can note that the r.h.s of Eq. (|?|) vanishes for < x < 
a since B(x) — in this range. Therefore, for x < a, one 
has the exact form P(x) = P(0) ex-p(—KP(0)x) which 
simply reflects the balance between the constant site cre- 
ation at x=0 and the depletion due to site recruitment 
in avalanches. Besides this simple result, an analytical 
determination of P(x) has eluded us. We have therefore 
solved numerically Eq. (|J) and (Q), taking advantage of 
the known form of P(x) for x < a. Given A(x), this 
determines the r.h.s of Eq. (Q) for a < x < 2a and thus 
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P(x) in this range. Continuing the process, the com- 
putation of P{x) consists of solving K linear inhomoge- 
neous differential equations. We have thus iterated the 
following process : solve Eq. M) for some P(0) and A, 
normalize P so that J P = 1, set P(0) to its new actual 
value, and use Eq. (||) to find the new A. The non-trivial 
satisfaction of relation (J3J) was used as a check of the 
computation. The computed P(x) for K = 4 is shown 
in Fig. |l] and agrees well with simulations results. For 
K = 2, similar agreement is obtained. The computed 
mean avalanche size agrees with the averaged avalanche 
size obtained from simulations, in the range of a where 
they can be compared, as shown in Fig. |^ for K = 2. For 
the largest a's, it was found necessary to use lattices of 
N = 2 10 4 and average over 3 10 5 avalanches to ensure 
that convergence to the steady state was reached and 
that the results were free from finite size effects. 
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FIG. 2. Mean avalanche size vs. a for K = 2 showing 
solutions of (^J) and (Q) with S = A + 1 (circles) and results 
of numerical simulations (filled square) . The line is a guide 
for the eyes. Insert : P{x) obtained by solving Eq. (|| fj]) for 
K = 2 and a = 0.45, 0.453, 0.456, 0.463. 

In both cases, the mean avalanche size grows very 
quickly when a approaches 1/K. In ref. H, similar 
results were interpreted as a divergence of the mean 
avalanche size at a c — 0.2255 for K — 4. With our nu- 
merical procedure, the solutions of Eq. (||) and ((?]) can, 
however, be reliably obtained up to a — 0.236, far above 
the purported critical a c . The corresponding distribu- 
tions P(x) are shown in Fig. ^ for K ~ 4 and do not dis- 
play any singularity at a = 0.2255. On the contrary, the 
four peaks of the distribution appear to sharpen smoothly 
as a increases. This leads us to think that they smoothly 
tend to 6 peaks as a — > 1/4. As shown in the insert of 
Fig. ||, a similar behavior is observed for K = 2. 

The presence of sharper and sharper peaks puts a 
heavy demand on numerical resolution and prevents a 
direct numerical approach of the limit a — ► 1/K (with 
our algorithm, at least). To analyze further this limit, 



we focus for simplicity on the case K — 2. For a = 1/2, 
P(x) is made of two 5 peaks located at x = and x = 1/2 
respectively while A{x) is simply a <5 peak at x=l. When 
a is close to 1/2, P(0) becomes large and the derivative 
term in Eq. (R) can balance the other terms only if P(x) 
has a fast variation on a scale 1/P(0). This is indeed the 
case of the exact form of P(x) near x = 0. We therefore 
search for A and P under the form, 



A(x)~ 



P(0) 2 



i[P(0)(a 



.1 l IWM.i J- — 277(0?) )1 

Fix) = P(0) exp(-2P(0)a;), x < " 
P(x) ~ i P(0) n[P(0)(x - a - 



(8) 



77(a))], x>a (9) 



where a and II are two functions to be determined which 
have been normalized so that their integrals equal one 
and which describe the broadening for a ^ 1/2 of the 8 
peaks at x = 1 and x = 1/2 respectively. The peak dis- 
placement 77(a) is supposed to tend to zero as a — > 1/2. 
Moreover, self-consistency requires that P(0)?7(a) —> 00 
when a tends toward 1/2. This allows to neglect the 
weight of B(x)/A near x = 1 (note that B(x) is the 
continuation of A(x) for x < 1) as supposed in (|9j). Sub- 
stituting the forms (||) and (||) in Eq. (||) and ([?]), we 
obtain at dominant order, 



a{x) 



+00 



U{x + 2C - u/2)a{u)du 



(10) 



/2x 
a(u)e u (x~u/2)du (11) 
-00 

where we have defined the constant C = lim Q ^ 1 / 2 (l/2 — 
a)P(0) and integrated the linear differential equation for 
II. Taking Fourier transforms of (|l0|) and (pT|), one ob- 
tains for a(uj) — J dxexp(iLux)a(x), the functional equa- 
tion 



d(w) 



exp(-2^C) „ 2 
(l-z^/2) 2 a( " /2) 



(12) 



We fix the translational symetry of Eq. ( |10| ) and ( |Tl| ) 
[a(x) — > a(x + a;o),II(a;) — > II(a: + Xo/2)] by impo sing 
J xa(x)dx = 0. Then, the unique solution of Eq. (|12| ) 
without a singularity at lo — is 



om = n 



exp(— 2iujC) 
\ (1 - iuj/2 n ) 2n 



(13) 



where convergence of the infinite product enforces C — 
1/2. One can check that the inverse Fourier transform 
a(x) of a(uj) is indeed a real function and is positive, as 
it should, since it is the convolution of the real positive 
functions 



v n (x) = 



2 n [2 n (x + l)f 
(2" - 1)! 



-2 n (x+l) , 



1) (14) 
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When x — > — oo, this allows to show that a(x) and H(x) 
tend extremely quickly towards zero, namely a(2x) ~ 
THx) ~ exp(— est 4~ x ). Comparing the two estimations 
(H) of near x = a, one obtains for the peak dis- 

placement P(0)r)(a) ~ | ln(l/2 - a)|/ln(4) as a -> 1/2. 
When x — > +00, the analytic expression ([l3|) (or di- 
rectly Eq. ( |i~0| ) and (pi])) shows that a(ic) and II(a;) tend 
toward zero as a;exp(— 2cc). Using this asymptotic be- 
havior to estimate P(l) gives that the mean avalanche 
size, P(0)/P(1), diverges like exp(cst/(l/2 — a)) when 
a — > 1/2. These predictions are compared in Fig. 4 to 
results obtained from the numerical solutions of Eq. (||) 
and (0) for K = 2 and different values of a. As shown in 
the insert, 1/P(0) vanishes linearly when a — ► 1/2 with a 
measured slope of 2.3 close to the analytical prediction of 
1/C = 2, the difference between the two being quite com- 
patible with higher order terms in 1/2 — a. The function 
a(x) obtained by taking the inverse Fourier transform of 
( |l3| ) compares well to rescaled plots of A(y) with 77(a) 
chosen so as to make the different maxima coincide [n0[ . 
Similar agreement is obtained between the analytic IU~x) 
and rescaled versions of P(y) around y — 1/2. Finally, 
the extremely rapid growth of the mean avalanche size 
when a — ► 1/2 agrees semi-quantitatively with the nu- 
merical results shown in Fig. 2 but is itself an obstacle to 
a precise numerical check of the predicted asymptotics. 
It makes it also difficult to avoid finite size effects in nu- 
merical simulations when a —> 1/ K given that the cut-off 
in the avalanche size distribution scales in a mean field 
manner as the square of the mean avalanche size || . 

In conclusion, evidence has been presented that the 
random neighbor OFC model has finite avalanches as 
soon as the model is non-conservative with a mean 
avalanche size which increases extremely quickly as the 
conservative limit is approached. It would be interest- 
ing to assess the generality of this phenomenon and its 
relevance for lattice models. 
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FIG. 4. a(x) (dashed line) and rescaled graphs A(y)a* /A* 
vs. x = (y — 1 — 2n(a))p (lines) for K = 2 and a = 0.45, 0.456 
and 0.463. a* /A* is the ratio of the curve maxima and 
p = A* /(a* J A); p/P(Q) = 1.145, 1.117, 1.089 for the graphs 
shown, approaching 1 as expected when a — + 1/2. Insert: 
1/P(0) versus (1/2 — a); the straight line fit has a slope of 2.3 



FIG. 3. Stationary distribution obtained by solving Eq. (|s| 
g) for K = 4 and a = 0.22, 0.224, 0.228, 0.232, 0.236. 
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